---
title: "Published peer review reports have higher informative
content than unpublished reports"
output: html_notebook
---

```{r}
library(dplyr)
library(plotly)
library(DescTools)
library(ggplot2)
library(dplyr)
library(gridExtra)
library(fs)
library(purrr)
library(readr)
library(reshape2)
library(PupillometryR)
library(stargazer)
library(modelsummary)
library(tidyr)
library(FactoMineR)
library(stringr) 
library(xtable)
library(gghalves)
library(glmmTMB)
library(lme4)
```

### Data loading
```{r}
dev_scores <- read.csv('./data_dev_scores_vf.csv')
gini_by_scientist <- read.csv('./data_gini_by_scientist_vf.csv')
gini_by_manuscript <- read.csv('./data_gini_by_manuscript_vf.csv')
readability_scores_ <- read_csv("./data_readability_scores_vf.csv")
```


### Figure 1: The distribution of the information score for published vs. unpublished peer review reports. 
```{r}
plot_data <- gini_by_manuscript %>% 
  mutate(series = "Gini Index") %>% 
  mutate(value = GiniIndex) %>% 
  mutate(`Peer Review Type` = as.character(PRType)) %>% 
  select(series, value, `Peer Review Type`) %>% 
  rbind(dev_scores %>% 
          mutate(series = "Info Score") %>% 
          mutate(value = Score) %>% 
          mutate(`Peer Review Type` = as.character(PRType)) %>% 
          select(series, value, `Peer Review Type`)) %>% 
  mutate(series = factor(series, levels = c("Gini Index", "Info Score" )))

violins_limits <- plot_data %>%
  group_by(series, `Peer Review Type`) %>%
  summarise(ymin = min(value), ymax = max(value))

medians <- plot_data %>%
  group_by(series, `Peer Review Type`) %>%
  summarise(median_value = median(value, na.rm = TRUE))

# Info Score
fig1 <- ggplot(data = plot_data %>% filter(series == "Info Score")) + 
  geom_flat_violin(aes(x = series, y = value, fill = `Peer Review Type`, color = `Peer Review Type`), 
                   position = "identity", alpha = 0.4, width = 1.8) +
  geom_segment(data = medians %>% filter(series == "Info Score"), 
               aes(x = as.numeric(series) - 1, xend = as.numeric(series) + 0.1,  # Ajuste de x
                   y = median_value, yend = median_value, color = `Peer Review Type`),
               linetype = "dashed", size = 0.8, show.legend = FALSE) +
  geom_text(data = medians %>% filter(series == "Info Score"), 
            aes(x = as.numeric(series), y = median_value, label = round(median_value, 2), 
                color = `Peer Review Type`, colour = "black"),
            position = position_dodge(width = 0.2), vjust = +28.5, hjust = -0.2, size = 6, show.legend = FALSE) +
  coord_flip() + 
  ylab("") +
  xlab("") + 
  scale_fill_grey() +
  scale_color_manual(values = c("Open" = "black", "Unpublished" = "darkgrey")) +
  theme(
    axis.text = element_text(size = 20),
    axis.text.x = element_text(margin = margin(t = -350)),
    legend.text = element_text(size = 20), 
    legend.title = element_text(size = 20),
    legend.position = "bottom",
    legend.margin = margin(t = -10),  
    panel.background = element_blank(),
    axis.title.x = element_blank()
  ) + 
  scale_x_discrete(labels = NULL) +
  guides(fill = guide_legend(reverse = TRUE), color = guide_legend(reverse = TRUE))  # Invertir el orden de la leyenda

fig1

```

### Figure 2: The distribution of the number of sentences per category for published vs. unpublished peer review reports.
```{r}
dims.names <- c("ExampleCount", "CriticismCount", "SuggestionAndSolutionCount", 
                "ImportanceAndRelevanceCount", "MaterialsAndMethodsCount", 
                "PraiseCount", "PresentationAndReportingCount", 
                "ResultsAndDiscussionCount")

data_long <- dev_scores %>%
  select(PRType, dims.names) %>%
  pivot_longer(cols = dims.names, 
               names_to = "count_type", 
               values_to = "count_value") %>%
  mutate(count_type = gsub("Count", "", count_type))       

# Reordenar los niveles de 'count_type'
data_long <- data_long %>%
  mutate(count_type = recode(count_type,
                             "SuggestionAndSolution" = "Suggestion And Solution",
                             "PresentationAndReporting" = "Presentation And Reporting",
                             "Praise" = "Praise",
                             "Example" = "Example",
                             "ImportanceAndRelevance" = "Importance And Relevance",
                             "ResultsAndDiscussion" = "Results And Discussion",
                             "MaterialsAndMethods" = "Materials And Methods",
                             "Criticism" = "Criticism"))

# Reordenar los niveles de 'count_type'
data_long <- data_long %>%
  mutate(count_type = factor(count_type, levels = c(
    "Suggestion And Solution",
    "Presentation And Reporting",
    "Praise",
    "Example",
    "Importance And Relevance",
    "Results And Discussion",
    "Materials And Methods",
    "Criticism"
  )))

fig2 <- ggplot(data_long, aes(x = count_value, y = count_type, fill = PRType)) + 
  geom_boxplot(position = position_dodge(width = 0.75),
               color = "black", size = 0.5) + 
  scale_fill_manual(values = c("Unpublished" = "lightgrey", "Open" = "darkgrey")) + 
  coord_cartesian(xlim = c(0, 13)) +  
  scale_x_continuous(breaks = seq(0, 13, by = 1)) + 
  theme_minimal() +
  labs(fill = "Peer Review Type",   x = "Number of Sentences", 
) +
  theme(
    axis.text.x = element_text(hjust = 1, size = 12),
    axis.text.y = element_text(size = 12),
    axis.title.x = element_text('number of sentences'),  
    axis.title.y = element_blank(),  
    legend.title = element_text(size = 14),
    legend.text = element_text(size = 12),
    legend.position = "bottom",
    plot.title = element_text(hjust = 0.5, size = 16),  
    panel.background = element_blank(),  
    plot.background = element_rect(fill = "white", color = NA),  
    panel.grid = element_blank(),  
    panel.border = element_blank()  
  ) +
  guides(fill = guide_legend(order = 1, reverse = TRUE, direction = "horizontal", nrow = 1))

fig2
```


### Figure 3: Distribution of the information score for signed (10,1673 reports) and unsigned (15,922 reports) published peer review reports.
```{r}
data_box <- dev_scores %>% 
  filter(PRType == 'Open') %>% 
  mutate(IsSigned = factor(IsSigned, levels = c("1", "0"))) %>% 
  mutate(key = 'IsSigned',
         series = IsSigned) %>% 
  mutate(key = factor(key, levels = c("IsSigned")))
  

fig3 <- ggplot(data_box, aes(x=series, y=Score)) + 
    geom_boxplot(width = 0.3) + 
    scale_fill_grey() +
    ylab("") +
    xlab("") + 
    theme(axis.text=element_text(size=22), legend.text=element_text(size=22), 
          legend.position="none", legend.title=element_blank(),
          strip.text.x = element_text(size = 22),
          axis.text.x = element_text(angle = 55, hjust = 1),
          panel.background = element_blank()) +
   scale_x_discrete(labels= c('Signed', 'Unsigned')) +
   labs(title = NULL)

fig3
```


### Figure 4: The distribution of the number of sentences for signed (10,1673 reports) and unsigned (15,922 reports) published peer review reports.
```{r}
data_box <- dev_scores %>% 
  filter(PRType == 'Open') %>% 
  mutate(IsSigned = factor(IsSigned, levels = c("1", "0"))) %>% 
  mutate(key = 'IsSigned',
         series = IsSigned) %>% 
  mutate(key = factor(key, levels = c("IsSigned"))) %>%
  select(PRType, dims.names, IsSigned) %>%
  pivot_longer(cols = dims.names, 
               names_to = "count_type", 
               values_to = "count_value") %>%
  mutate(count_type = gsub("Count", "", count_type)) %>%
  mutate(count_type = recode(count_type,
                             "SuggestionAndSolution" = "Suggestion And Solution",
                             "PresentationAndReporting" = "Presentation And Reporting",
                             "Praise" = "Praise",
                             "Example" = "Example",
                             "ImportanceAndRelevance" = "Importance And Relevance",
                             "ResultsAndDiscussion" = "Results And Discussion",
                             "MaterialsAndMethods" = "Materials And Methods",
                             "Criticism" = "Criticism"),
         count_type = factor(count_type, levels = c(
           "Suggestion And Solution",
           "Presentation And Reporting",
           "Praise",
           "Example",
           "Importance And Relevance",
           "Results And Discussion",
           "Materials And Methods",
           "Criticism"
         )),
         IsSigned = recode(IsSigned, "1" = "Signed", "0" = "Unsigned"))  # recodifica aquí

median_values <- data_box %>%
  filter(IsSigned == "Unsigned") %>%
  group_by(count_type) %>%
  summarise(median_value = median(count_value, na.rm = TRUE), .groups = "drop") %>%
  arrange(desc(median_value))

data_box$count_type <- factor(data_box$count_type, levels = median_values$count_type)

fig4 <- ggplot(data_box, aes(x = count_value, y = count_type, fill = IsSigned)) + 
  geom_boxplot(position = position_dodge(width = 0.75),
               color = "black", size = 0.5) + 
  scale_fill_manual(values = c("Signed" = "lightgrey", "Unsigned" = "darkgrey")) +
  coord_cartesian(xlim = c(0, 13)) +  
  scale_x_continuous(breaks = seq(0, 13, by = 1)) + 
  theme_minimal() +
  labs(fill = "Reviewer identity", x = "Number of Sentences") +
  theme(
    axis.text.x = element_text(hjust = 1, size = 12),
    axis.text.y = element_text(size = 12),
    axis.title.x = element_text('Number of Sentences'),  
    axis.title.y = element_blank(),  
    legend.title = element_text(size = 14),
    legend.text = element_text(size = 12),
    legend.position = "bottom",
    plot.title = element_text(hjust = 0.5, size = 16),  
    panel.background = element_blank(),  
    plot.background = element_rect(fill = "white", color = NA),  
    panel.grid = element_blank(),  
    panel.border = element_blank()  
  ) +
  guides(fill = guide_legend(order = 1, reverse = TRUE, direction = "horizontal", nrow = 1))


fig4
```


### Figure 5: The dispersion of the information score (Gini index) for published vs. unpublished peer review reports. 
```{r}
# Gini Index
fig_gini_index <- ggplot(data = plot_data %>% filter(series == "Gini Index")) + 
  geom_flat_violin(aes(x = series, y = value, fill = `Peer Review Type`, color = `Peer Review Type`), 
                   position = "identity", alpha = 0.4, width = 1.8) +
  geom_segment(data = medians %>% filter(series == "Gini Index"), 
               aes(x = as.numeric(series), xend = as.numeric(series) + 1, 
                   y = median_value, yend = median_value, color = `Peer Review Type`),
               linetype = "dashed", size = 0.8, show.legend = FALSE) +
  geom_text(data = medians %>% filter(series == "Gini Index"), 
            aes(x = series, y = median_value, label = round(median_value, 2), 
                color = `Peer Review Type`, colour="black"),
            position = position_dodge(width = 0.2), vjust = -3.5, hjust= -0.2, size = 6, show.legend = FALSE) +
  coord_flip() + 
  ylab("") +
  xlab("") + 
  scale_fill_grey() +
  scale_color_manual(values = c("Open" = "black", "Unpublished" = "darkgrey")) +
  theme(
    axis.text = element_text(size = 20),
    axis.text.x = element_text(margin = margin(t = -350)),
    legend.text = element_text(size = 20), 
    legend.title = element_text(size = 20),
    legend.position = "bottom",
    legend.margin = margin(t = -10),  
    panel.background = element_blank(),
    axis.title.x = element_blank()
  ) + scale_x_discrete(labels = NULL)


```

```{r}
# CI
ci_ginimanuscript <- t.test(gini_by_manuscript[gini_by_manuscript$PRType == 'Unpublished',]$GiniIndex, 
                gini_by_manuscript[gini_by_manuscript$PRType != 'Unpublished',]$GiniIndex)$conf.int

ci_score <- t.test(dev_scores[dev_scores$PRType == 'Unpublished',]$Score, 
                dev_scores[dev_scores$PRType != 'Unpublished',]$Score)$conf.int

ci_ex <- t.test(dev_scores[dev_scores$PRType == 'Unpublished',]$Example, 
                dev_scores[dev_scores$PRType != 'Unpublished',]$Example)$conf.int

ci_cri <- t.test(dev_scores[dev_scores$PRType == 'Unpublished',]$Criticism, 
                 dev_scores[dev_scores$PRType != 'Unpublished',]$Criticism)$conf.int

ci_mat <- t.test(dev_scores[dev_scores$PRType == 'Unpublished',]$MaterialsAndMethods, 
                 dev_scores[dev_scores$PRType != 'Unpublished',]$MaterialsAndMethods)$conf.int

ci_su <- t.test(dev_scores[dev_scores$PRType == 'Unpublished',]$SuggestionAndSolution, 
                dev_scores[dev_scores$PRType != 'Unpublished',]$SuggestionAndSolution)$conf.int

ci_imp <- t.test(dev_scores[dev_scores$PRType == 'Unpublished',]$ImportanceAndRelevance, 
                 dev_scores[dev_scores$PRType != 'Unpublished',]$ImportanceAndRelevance)$conf.int

ci_pra <- t.test(dev_scores[dev_scores$PRType == 'Unpublished',]$Praise, 
                 dev_scores[dev_scores$PRType != 'Unpublished',]$Praise)$conf.int

ci_pre <- t.test(dev_scores[dev_scores$PRType == 'Unpublished',]$PresentationAndReporting, 
                 dev_scores[dev_scores$PRType != 'Unpublished',]$PresentationAndReporting)$conf.int

ci_re <- t.test(dev_scores[dev_scores$PRType == 'Unpublished',]$ResultsAndDiscussion, 
                dev_scores[dev_scores$PRType != 'Unpublished',]$ResultsAndDiscussion)$conf.int


# quantiles 25% y 75%
iqr_values <- list(
  GiniManuscript = list(
    Unpublished = quantile(gini_by_manuscript[gini_by_manuscript$PRType == 'Unpublished',]$GiniIndex, probs = c(0.25, 0.75)),
    Open = quantile(gini_by_manuscript[gini_by_manuscript$PRType != 'Unpublished',]$GiniIndex, probs = c(0.25, 0.75))
  ),
  Score = list(
    Unpublished = quantile(dev_scores[dev_scores$PRType == 'Unpublished',]$Score, probs = c(0.25, 0.75)),
    Open = quantile(dev_scores[dev_scores$PRType != 'Unpublished',]$Score, probs = c(0.25, 0.75))
  ),
  NumberSentences = list(
    Unpublished = quantile(gini_by_manuscript[gini_by_manuscript$PRType == 'Unpublished',]$GiniIndex, probs = c(0.25, 0.75)),
    Open = quantile(gini_by_manuscript[gini_by_manuscript$PRType != 'Unpublished',]$GiniIndex, probs = c(0.25, 0.75))
  ),
  Example = list(
    Unpublished = quantile(dev_scores[dev_scores$PRType == 'Unpublished',]$Example, probs = c(0.25, 0.75)),
    Open = quantile(dev_scores[dev_scores$PRType != 'Unpublished',]$Example, probs = c(0.25, 0.75))
  ),
  Criticism = list(
    Unpublished = quantile(dev_scores[dev_scores$PRType == 'Unpublished',]$Criticism, probs = c(0.25, 0.75)),
    Open = quantile(dev_scores[dev_scores$PRType != 'Unpublished',]$Criticism, probs = c(0.25, 0.75))
  ),
  MaterialsAndMethods = list(
    Unpublished = quantile(dev_scores[dev_scores$PRType == 'Unpublished',]$MaterialsAndMethods, probs = c(0.25, 0.75)),
    Open = quantile(dev_scores[dev_scores$PRType != 'Unpublished',]$MaterialsAndMethods, probs = c(0.25, 0.75))
  ),
  SuggestionAndSolution = list(
    Unpublished = quantile(dev_scores[dev_scores$PRType == 'Unpublished',]$SuggestionAndSolution, probs = c(0.25, 0.75)),
    Open = quantile(dev_scores[dev_scores$PRType != 'Unpublished',]$SuggestionAndSolution, probs = c(0.25, 0.75))
  ),
  ImportanceAndRelevance = list(
    Unpublished = quantile(dev_scores[dev_scores$PRType == 'Unpublished',]$ImportanceAndRelevance, probs = c(0.25, 0.75)),
    Open = quantile(dev_scores[dev_scores$PRType != 'Unpublished',]$ImportanceAndRelevance, probs = c(0.25, 0.75))
  ),
  Praise = list(
    Unpublished = quantile(dev_scores[dev_scores$PRType == 'Unpublished',]$Praise, probs = c(0.25, 0.75)),
    Open = quantile(dev_scores[dev_scores$PRType != 'Unpublished',]$Praise, probs = c(0.25, 0.75))
  ),
  PresentationAndReporting = list(
    Unpublished = quantile(dev_scores[dev_scores$PRType == 'Unpublished',]$PresentationAndReporting, probs = c(0.25, 0.75)),
    Open = quantile(dev_scores[dev_scores$PRType != 'Unpublished',]$PresentationAndReporting, probs = c(0.25, 0.75))
  ),
  ResultsAndDiscussion = list(
    Unpublished = quantile(dev_scores[dev_scores$PRType == 'Unpublished',]$ResultsAndDiscussion, probs = c(0.25, 0.75)),
    Open = quantile(dev_scores[dev_scores$PRType != 'Unpublished',]$ResultsAndDiscussion, probs = c(0.25, 0.75))
  )
)
```

### Table 1: Table 1: Generalized linear mixed model on the information score by type of peer review, journal’s quartile of impact factor, reviewer gender and region of affiliation with journal-level random effects
```{r}
dev_scores$JournalID <- as.factor(dev_scores$JournalID)

dev_scores_ <- dev_scores %>%
  mutate(GenderBin = if_else(Gender == "Missing Gender", NA, Gender)) 

dev_scores_$EuNaAu <- as.factor(dev_scores_$EuNaAu)
dev_scores_$EuNaAu <- relevel(dev_scores_$EuNaAu, ref = "TRUE")


control_glmer <- glmerControl(
  optimizer = "bobyqa",             
  optCtrl = list(maxfun = 1e5)      
)

model_mixed <- glmer(
  Score ~ PRType + IFQuartile + GenderBin * EuNaAu + (1 | JournalID),
  data = dev_scores_,
  family = Gamma(link = "log"), 
  control = control_glmer
)

stargazer(model_mixed, single.row=TRUE)

modelsummary(model_mixed, 
             stars = TRUE,
             exponentiate = FALSE,
             estimate = "{estimate} [{conf.low}, {conf.high}] ({std.error})", 
             statistic = "p={p.value}{stars}")
```

### Table 2: Generalized linear model of the dispersion of the information score (Gini index) by type of peer review, journal’s quartile of impact factor, reviewer gender and region of affiliation
```{r}
gini_by_manuscript$EuNaAuAgg <- as.factor(gini_by_manuscript$EuNaAuAgg)
gini_by_manuscript$EuNaAuAgg <- relevel(gini_by_manuscript$EuNaAuAgg, ref = "All EuNaAu")

table3 = glm(GiniIndex ~ PRType+IFQuartile+GenderAgg*EuNaAuAgg , 
             data=gini_by_manuscript)
stargazer(table3, single.row=TRUE)

modelsummary(table3, 
             exponentiate = FALSE,
             stars = TRUE,
             estimate="{estimate} [{conf.low}, {conf.high}] ({std.error})", 
             statistic = "p={p.value}{stars}")

```


### Table 3: Generalized linear mixed model on the report readibility with standardised (z-scored) coefficients and percentage of change by type of peer review, journal impact factor quartile, reviewer gender and region of affiliation with journal-level random effects
```{r}
readability_scores <- readability_scores_ %>% 
  group_by(ManuscriptID) %>% 
  mutate(nb=n(), numberManuscripts=n()) %>% 
  mutate(GenderText = paste0(Gender, collapse = "-")) %>% 
  mutate(EuNaAuText = paste0(ContinentalRegion, collapse = "-")) %>% 
  mutate(
    GenderAgg = case_when(
      all(Gender == "Women") ~ "All women",
      all(Gender == "Men") ~ "All men",
      TRUE ~ "Mixed"
    ),
    EuNaAuAgg = case_when(
      all(ContinentalRegion %in% c("Europe", "North America", "Oceania")) ~ "All EuNaAu",
      all(!ContinentalRegion %in% c("Europe", "North America", "Oceania")) ~ "All rest",
      TRUE ~ "Mixed"
    )
  )  %>% 
  ungroup() %>% 
  mutate(
    EuNaAu = ifelse(ContinentalRegion %in% c("Europe", "North America", "Oceania"), TRUE, FALSE),
    flesch_interpret = case_when(
      flesch_ease >= 90 ~ "very easy - 5th grade",
      flesch_ease >= 80 ~ "easy - 5th grade",
      flesch_ease >= 70 ~ "fairly easy - 5th grade",
      flesch_ease >= 60 ~ "standard - 8th-9th grade",
      flesch_ease >= 50 ~ "fairly difficult - 10th-12th grade",
      flesch_ease >= 30 ~ "difficult - college",
      TRUE ~ "very difficult - college graduate"),
    smog_interpret = case_when(
      smog == 6 ~ "Very easy to read – 6th grade",
      smog == 7 ~ "Easy to read – 7th grade",
      smog == 8 ~ "Fairly easy to read – 8th grade",
      smog == 9 ~ "Plain English – 9th grade",
      smog == 10 ~ "Fairly difficult to read – 10th grade",
      smog == 11 ~ "Difficult to read – 11th grade",
      smog == 12 ~ "Very difficult to read – 12th grade",
      smog >= 13 ~ "College level and above – College and beyond",
      TRUE ~ "College level and above – College and beyond"
    )
  ) 
  
readability_scores2 <- readability_scores %>% 
  filter(ManuscriptID %in% dev_scores$ManuscriptID)

readability_scores$PRType <- factor(readability_scores$PRType, levels = c("Unpublished", "Open"))
readability_scores$IFQuartile <- factor(readability_scores$IFQuartile, levels = c("Q1","Q2","Q3","Q4+NI"))
readability_scores$EuNaAu <- factor(readability_scores$EuNaAu, levels = c("TRUE", "FALSE"))

readability_scores <- readability_scores %>% 
  mutate(GenderBin = if_else(Gender == "Missing Gender", NA, Gender)) 

readability_scores$GenderBin <- factor(readability_scores$GenderBin)


control_lmer <- lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e5))

metrics <- c("hf_readability_grade_mdeberta")

results <- lapply(metrics, function(metric) {
  
  df <- readability_scores %>%
    mutate(metric_z = scale(.data[[metric]]))
  
  model <- lmer(
    metric_z ~ PRType + IFQuartile + GenderBin * EuNaAu + (1 | Journal),
    data = df,
    control = control_lmer
  )
  
  stargazer(model, single.row=TRUE)

  mu <- mean(df[[metric]], na.rm = TRUE)
  sigma <- sd(df[[metric]], na.rm = TRUE)
  
  
  tidy(model, effects = "fixed") %>%
    select(term, estimate, std.error) %>%
    mutate(
      metric = metric,
      delta_metric = estimate * sigma,
      percent_change = (delta_metric / mu) * 100
    )
})

coef_all_metrics <- bind_rows(results)
coef_all_metrics

```

